<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"
                "http://www.w3.org/TR/REC-html40/loose.dtd">
<html>
<head>
  <title>Description of demgmm4</title>
  <meta name="keywords" content="demgmm4">
  <meta name="description" content="DEMGMM4 Demonstrate density modelling with a Gaussian mixture model.">
  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">
  <meta name="generator" content="m2html &copy; 2003 Guillaume Flandin">
  <meta name="robots" content="index, follow">
  <link type="text/css" rel="stylesheet" href="../../m2html.css">
</head>
<body>
<a name="_top"></a>
<div><a href="../../menu.html">Home</a> &gt;  <a href="#">ReBEL-0.2.7</a> &gt; <a href="#">netlab</a> &gt; demgmm4.m</div>

<!--<table width="100%"><tr><td align="left"><a href="../../menu.html"><img alt="<" border="0" src="../../left.png">&nbsp;Master index</a></td>
<td align="right"><a href="menu.html">Index for .\ReBEL-0.2.7\netlab&nbsp;<img alt=">" border="0" src="../../right.png"></a></td></tr></table>-->

<h1>demgmm4
</h1>

<h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>DEMGMM4 Demonstrate density modelling with a Gaussian mixture model.</strong></div>

<h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="box"><strong>This is a script file. </strong></div>

<h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre class="comment">DEMGMM4 Demonstrate density modelling with a Gaussian mixture model.

    Description
     The problem consists of modelling data generated by a mixture of
    three Gaussians in 2 dimensions with a mixture model using full
    covariance matrices.  The priors are 0.3, 0.5 and 0.2; the centres
    are (2, 3.5), (0, 0) and (0,2); the variances are (0.16, 0.64) axis
    aligned, (0.25, 1) rotated by 30 degrees and the identity matrix. The
    first figure contains a scatter plot of the data.

    A Gaussian mixture model with three components is trained using EM.
    The parameter vector is printed before training and after training.
    The user should press any key to continue at these points.  The
    parameter vector consists of priors (the column), and centres (given
    as (x, y) pairs as the next two columns).  The covariance matrices
    are printed separately.

    The second figure is a 3 dimensional view of the density function,
    while the third shows the axes of the 1-standard deviation ellipses
    for the three components of the mixture model.

    See also
    <a href="gmm.html" class="code" title="function mix = gmm(dim, ncentres, covar_type, ppca_dim)">GMM</a>, <a href="gmminit.html" class="code" title="function mix = gmminit(mix, x, options)">GMMINIT</a>, <a href="gmmem.html" class="code" title="function [mix, options, errlog] = gmmem(mix, x, options)">GMMEM</a>, <a href="gmmprob.html" class="code" title="function prob = gmmprob(mix, x)">GMMPROB</a>, <a href="gmmunpak.html" class="code" title="function mix = gmmunpak(mix, p)">GMMUNPAK</a></pre></div>

<!-- crossreference -->
<h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
This function calls:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="gmm.html" class="code" title="function mix = gmm(dim, ncentres, covar_type, ppca_dim)">gmm</a>	GMM	Creates a Gaussian mixture model with specified architecture.</li><li><a href="gmmem.html" class="code" title="function [mix, options, errlog] = gmmem(mix, x, options)">gmmem</a>	GMMEM	EM algorithm for Gaussian mixture model.</li><li><a href="gmminit.html" class="code" title="function mix = gmminit(mix, x, options)">gmminit</a>	GMMINIT Initialises Gaussian mixture model from data</li><li><a href="gmmprob.html" class="code" title="function prob = gmmprob(mix, x)">gmmprob</a>	GMMPROB Computes the data probability for a Gaussian mixture model.</li></ul>
This function is called by:
<ul style="list-style-image:url(../../matlabicon.gif)">
<li><a href="demnlab.html" class="code" title="function demnlab(action);">demnlab</a>	DEMNLAB A front-end Graphical User Interface to the demos</li></ul>
<!-- crossreference -->


<h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>
<div class="fragment"><pre>0001 <span class="comment">%DEMGMM4 Demonstrate density modelling with a Gaussian mixture model.</span>
0002 <span class="comment">%</span>
0003 <span class="comment">%    Description</span>
0004 <span class="comment">%     The problem consists of modelling data generated by a mixture of</span>
0005 <span class="comment">%    three Gaussians in 2 dimensions with a mixture model using full</span>
0006 <span class="comment">%    covariance matrices.  The priors are 0.3, 0.5 and 0.2; the centres</span>
0007 <span class="comment">%    are (2, 3.5), (0, 0) and (0,2); the variances are (0.16, 0.64) axis</span>
0008 <span class="comment">%    aligned, (0.25, 1) rotated by 30 degrees and the identity matrix. The</span>
0009 <span class="comment">%    first figure contains a scatter plot of the data.</span>
0010 <span class="comment">%</span>
0011 <span class="comment">%    A Gaussian mixture model with three components is trained using EM.</span>
0012 <span class="comment">%    The parameter vector is printed before training and after training.</span>
0013 <span class="comment">%    The user should press any key to continue at these points.  The</span>
0014 <span class="comment">%    parameter vector consists of priors (the column), and centres (given</span>
0015 <span class="comment">%    as (x, y) pairs as the next two columns).  The covariance matrices</span>
0016 <span class="comment">%    are printed separately.</span>
0017 <span class="comment">%</span>
0018 <span class="comment">%    The second figure is a 3 dimensional view of the density function,</span>
0019 <span class="comment">%    while the third shows the axes of the 1-standard deviation ellipses</span>
0020 <span class="comment">%    for the three components of the mixture model.</span>
0021 <span class="comment">%</span>
0022 <span class="comment">%    See also</span>
0023 <span class="comment">%    GMM, GMMINIT, GMMEM, GMMPROB, GMMUNPAK</span>
0024 <span class="comment">%</span>
0025 
0026 <span class="comment">%    Copyright (c) Ian T Nabney (1996-2001)</span>
0027 
0028 
0029 <span class="comment">% Generate the data</span>
0030 
0031 ndata = 500;
0032 
0033 <span class="comment">% Fix the seeds for reproducible results</span>
0034 randn(<span class="string">'state'</span>, 42);
0035 rand(<span class="string">'state'</span>, 42);
0036 data = randn(ndata, 2);
0037 prior = [0.3 0.5 0.2];
0038 <span class="comment">% Mixture model swaps clusters 1 and 3</span>
0039 datap = [0.2 0.5 0.3];
0040 datac = [0 2; 0 0; 2 3.5];
0041 datacov = repmat(eye(2), [1 1 3]);
0042 data1 = data(1:prior(1)*ndata,:);
0043 data2 = data(prior(1)*ndata+1:(prior(2)+prior(1))*ndata, :);
0044 data3 = data((prior(1)+prior(2))*ndata +1:ndata, :);
0045 
0046 <span class="comment">% First cluster has axis aligned variance and centre (2, 3.5)</span>
0047 data1(:, 1) = data1(:, 1)*0.4 + 2.0;
0048 data1(:, 2) = data1(:, 2)*0.8 + 3.5;
0049 datacov(:, :, 3) = [0.4*0.4 0; 0 0.8*0.8];
0050 
0051 <span class="comment">% Second cluster has variance axes rotated by 30 degrees and centre (0, 0)</span>
0052 rotn = [cos(pi/6) -sin(pi/6); sin(pi/6) cos(pi/6)];
0053 data2(:,1) = data2(:, 1)*0.5;
0054 data2 = data2*rotn;
0055 datacov(:, :, 2) = rotn' * [0.25 0; 0 1] * rotn;
0056 
0057 <span class="comment">% Third cluster is at (0,2)</span>
0058 data3 = data3 + repmat([0 2], prior(3)*ndata, 1);
0059 
0060 <span class="comment">% Put the dataset together again</span>
0061 data = [data1; data2; data3];
0062 
0063 clc
0064 disp(<span class="string">'This demonstration illustrates the use of a Gaussian mixture model'</span>)
0065 disp(<span class="string">'with full covariance matrices to approximate the unconditional '</span>)
0066 disp(<span class="string">'probability density of data in a two-dimensional space.'</span>)
0067 disp(<span class="string">'We begin by generating the data from a mixture of three Gaussians and'</span>)
0068 disp(<span class="string">'plotting it.'</span>)
0069 disp(<span class="string">' '</span>)
0070 disp(<span class="string">'The first cluster has axis aligned variance and centre (0, 2).'</span>)
0071 disp(<span class="string">'The second cluster has variance axes rotated by 30 degrees'</span>)
0072 disp(<span class="string">'and centre (0, 0).  The third cluster has unit variance and centre'</span>)
0073 disp(<span class="string">'(2, 3.5).'</span>)
0074 disp(<span class="string">' '</span>)
0075 disp(<span class="string">'Press any key to continue.'</span>)
0076 pause
0077 
0078 fh1 = figure;
0079 plot(data(:, 1), data(:, 2), <span class="string">'o'</span>)
0080 set(gca, <span class="string">'Box'</span>, <span class="string">'on'</span>)
0081 
0082 <span class="comment">% Set up mixture model</span>
0083 ncentres = 3;
0084 input_dim = 2;
0085 mix = <a href="gmm.html" class="code" title="function mix = gmm(dim, ncentres, covar_type, ppca_dim)">gmm</a>(input_dim, ncentres, <span class="string">'full'</span>);
0086 
0087 <span class="comment">% Initialise the model parameters from the data</span>
0088 options = foptions;
0089 options(14) = 5;    <span class="comment">% Just use 5 iterations of k-means in initialisation</span>
0090 mix = <a href="gmminit.html" class="code" title="function mix = gmminit(mix, x, options)">gmminit</a>(mix, data, options);
0091 
0092 <span class="comment">% Print out model</span>
0093 clc
0094 disp(<span class="string">'The mixture model has three components and full covariance'</span>)
0095 disp(<span class="string">'matrices.  The model parameters after initialisation using the'</span>)
0096 disp(<span class="string">'k-means algorithm are as follows'</span>)
0097 disp(<span class="string">'    Priors        Centres'</span>)
0098 disp([mix.priors' mix.centres])
0099 disp(<span class="string">'Covariance matrices are'</span>)
0100 disp(mix.covars)
0101 disp(<span class="string">'Press any key to continue.'</span>)
0102 pause
0103 
0104 <span class="comment">% Set up vector of options for EM trainer</span>
0105 options = zeros(1, 18);
0106 options(1)  = 1;        <span class="comment">% Prints out error values.</span>
0107 options(14) = 50;        <span class="comment">% Number of iterations.</span>
0108 
0109 disp(<span class="string">'We now train the model using the EM algorithm for 50 iterations.'</span>)
0110 disp(<span class="string">' '</span>)
0111 disp(<span class="string">'Press any key to continue.'</span>)
0112 pause
0113 [mix, options, errlog] = <a href="gmmem.html" class="code" title="function [mix, options, errlog] = gmmem(mix, x, options)">gmmem</a>(mix, data, options);
0114 
0115 <span class="comment">% Print out model</span>
0116 disp(<span class="string">' '</span>)
0117 disp(<span class="string">'The trained model has priors and centres:'</span>)
0118 disp(<span class="string">'    Priors        Centres'</span>)
0119 disp([mix.priors' mix.centres])
0120 disp(<span class="string">'The data generator has priors and centres'</span>)
0121 disp(<span class="string">'    Priors        Centres'</span>)
0122 disp([datap' datac])
0123 disp(<span class="string">'Model covariance matrices are'</span>)
0124 disp(mix.covars(:, :, 1))
0125 disp(mix.covars(:, :, 2))
0126 disp(mix.covars(:, :, 3))
0127 disp(<span class="string">'Data generator covariance matrices are'</span>)
0128 disp(datacov(:, :, 1))
0129 disp(datacov(:, :, 2))
0130 disp(datacov(:, :, 3))
0131 disp(<span class="string">'Note the close correspondence between these parameters and those'</span>)
0132 disp(<span class="string">'of the distribution used to generate the data.  The match for'</span>)
0133 disp(<span class="string">'covariance matrices is not that close, but would be improved with'</span>)
0134 disp(<span class="string">'more iterations of the training algorithm.'</span>)
0135 disp(<span class="string">' '</span>)
0136 disp(<span class="string">'Press any key to continue.'</span>)
0137 pause
0138 
0139 clc
0140 disp(<span class="string">'We now plot the density given by the mixture model as a surface plot.'</span>)
0141 disp(<span class="string">' '</span>)
0142 disp(<span class="string">'Press any key to continue.'</span>)
0143 pause
0144 
0145 <span class="comment">% Plot the result</span>
0146 x = -4.0:0.2:5.0;
0147 y = -4.0:0.2:5.0;
0148 [X, Y] = meshgrid(x,y);
0149 X = X(:);
0150 Y = Y(:);
0151 grid = [X Y];
0152 Z = <a href="gmmprob.html" class="code" title="function prob = gmmprob(mix, x)">gmmprob</a>(mix, grid);
0153 Z = reshape(Z, length(x), length(y));
0154 c = mesh(x, y, Z);
0155 hold on
0156 title(<span class="string">'Surface plot of probability density'</span>)
0157 hold off
0158 drawnow
0159 
0160 clc
0161 disp(<span class="string">'The final plot shows the centres and widths, given by one standard'</span>)
0162 disp(<span class="string">'deviation, of the three components of the mixture model.  The axes'</span>)
0163 disp(<span class="string">'of the ellipses of constant density are shown.'</span>)
0164 disp(<span class="string">' '</span>)
0165 disp(<span class="string">'Press any key to continue.'</span>)
0166 pause
0167 
0168 <span class="comment">% Try to calculate a sensible position for the second figure, below the first</span>
0169 fig1_pos = get(fh1, <span class="string">'Position'</span>);
0170 fig2_pos = fig1_pos;
0171 fig2_pos(2) = fig2_pos(2) - fig1_pos(4) - 30;
0172 fh2 = figure(<span class="string">'Position'</span>, fig2_pos);
0173 
0174 h3 = plot(data(:, 1), data(:, 2), <span class="string">'bo'</span>);
0175 axis equal;
0176 hold on
0177 title(<span class="string">'Plot of data and covariances'</span>)
0178 <span class="keyword">for</span> i = 1:ncentres
0179   [v,d] = eig(mix.covars(:,:,i));
0180   <span class="keyword">for</span> j = 1:2
0181     <span class="comment">% Ensure that eigenvector has unit length</span>
0182     v(:,j) = v(:,j)/norm(v(:,j));
0183     start=mix.centres(i,:)-sqrt(d(j,j))*(v(:,j)');
0184     endpt=mix.centres(i,:)+sqrt(d(j,j))*(v(:,j)');
0185     linex = [start(1) endpt(1)];
0186     liney = [start(2) endpt(2)];
0187     line(linex, liney, <span class="string">'Color'</span>, <span class="string">'k'</span>, <span class="string">'LineWidth'</span>, 3)
0188   <span class="keyword">end</span>
0189   <span class="comment">% Plot ellipses of one standard deviation</span>
0190   theta = 0:0.02:2*pi;
0191   x = sqrt(d(1,1))*cos(theta);
0192   y = sqrt(d(2,2))*sin(theta);
0193   <span class="comment">% Rotate ellipse axes</span>
0194   ellipse = (v*([x; y]))';
0195   <span class="comment">% Adjust centre</span>
0196   ellipse = ellipse + ones(length(theta), 1)*mix.centres(i,:);
0197   plot(ellipse(:,1), ellipse(:,2), <span class="string">'r-'</span>);
0198 <span class="keyword">end</span>
0199 hold off
0200 
0201 disp(<span class="string">'Note how the data cluster positions and widths are captured by'</span>)
0202 disp(<span class="string">'the mixture model.'</span>)
0203 disp(<span class="string">' '</span>)
0204 disp(<span class="string">'Press any key to end.'</span>)
0205 pause
0206 
0207 close(fh1);
0208 close(fh2);
0209 clear all; 
0210</pre></div>
<hr><address>Generated on Tue 26-Sep-2006 10:36:21 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/">m2html</a></strong> &copy; 2003</address>
</body>
</html>